Robust gene regulation: 
Deterministic dynamics from asynchronous networks with delay 
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We compare asynchronous vs. synchronous update of discrete dynamical networks and find that 
a simple time delay in the nodes may induce a reproducible deterministic dynamics even in the case 
of asynchronous update in random order. In particular we observe that the dynamics under syn- 
chronous parallel update can be reproduced accurately under random asynchronous serial update for 
a large class of networks. This mechanism points at a possible general principle of how computation 
in gene regulation networks can be kept in a quasi-deterministic "clockwork mode" in spite of the 
absence of a central clock. A delay similar to the one occurring in gene regulation causes synchro- 
nization in the model. Stability under asynchronous dynamics disfavors topologies containing loops, 
comparing well with the observed strong suppression of loops in biological regulatory networks. 

PACS numbers: 87.16.Yc,05.45.Xt,89.75.Hc,05.65.+b 
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Erwin Schrodinger in his lecture "What is life?" held 
in 1943 was one of the first to notice that the in- 
formation processing performed in the living cell has 
to be extremely robust and therefore requires a quasi- 
deterministic dynamics (which he called "clockwork 
mode"). The discovery of a "digital" storage medium for 
the genetic information, the double-stranded DNA, con- 
firmed one important part of this picture. Today, new 
experimental techniques allow to observe the dynamics 
of regulatory genes in great detail, which motivates us 
to reconsider the other, dynamical part of Schrodinger 's 
picture of a "clockwork mode" . While the dynamical ele- 
ments of gene regulation often are known in great detail, 
the complex dynamical patterns of the vast network of 
interacting regulatory genes, while highly reproducible 
between identical cells and organisms under similar con- 
ditions, are largely not understood. Most remarkably, 
these virtually deterministic activation patterns are often 
generated by asynchronous genetic switches without any 
central clock. In this Letter we address this astonishing 
fact with a toy model of gene regulation and study the 
conditions of when deterministic dynamics could occur in 
asynchronous circuits. Let us start from the observed dy- 
namics of small circuits of regulatory genes, then derive 
a discrete dynamical model gene, followed by a study of 
networks of such genetic switches, with a focus on com- 
paring their asynchronous and synchronous dynamics. 

Recently, several small gene regulation circuits have 
been described in terms of a detailed picture of their dy- 
namics 0, 0, 0, H, El A particularly simple motif is the 
single, self-regulating gene 0,0 that allows for a detailed 
modeling of its dynamics. A set of two differential equa- 
tions, for the temporal evolution of the concentrations of 
messenger RNA and protein, respectively, and an explicit 
time delay for transmission delay provide a quantitative 
model for the observed dynamics in this minimal circuit 



The equations of this model take the basic form 
dc 



dt 
db 

d* 



= a[f(s(t - *)) - c(t)] 

= pw) - m 



(i) 

(2) 



for the the dynamics of the concentrations c of mRNA 
and b of protein, with some non-linear transmission func- 
tion /(s) of an input signal s, a time delay and the 
time constants a and f3. In order to define a minimal 
discrete gene model let us keep the basic features (delay, 
low pass filter characteristics), omit the second filter, and 
write the difference equation for one gene i as 



Acj = a[f(si(t - -&)) - Cj(i)]At 



(3) 



The non-linear function / is typically a steep sigmoid. 
We approximate it as a step function O with 8(s) = 
for s < and O(s) = 1 otherwise. Rescaling time with 
e = a At and r = $/ At this reads 



Ac, = e[Q( Si (t - t)) - d(t)] 



(4) 



For simplicity let us update c$ by equidistant steps ac- 
cording to 

{+e, if Si(t — t) > and Cj < 1 — e 
— e, if Si(t — t) < and Cj > e (5) 
0, otherwise 

The coupling between nodes is defined by 

s i(t) = Wi X A t ) ~ a i I ( 6 ) 

3 

with discrete output states Xj (t) of the nodes defined as 
Xj (t) = Q( Cj (t)-l/2) . (7) 
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The influence of node j on node i can be activating 
(wij = 1), inhibitory = —1), or absent (my = 0). A 
constant bias cij is assigned to each node. 

In the following let us consider a network model ol 
such nodes. Consider N nodes with concentration vari- 
ables Ci, state variables Xj, biases a, and a coupling ma- 
trix (wij). Given initial values x,(0) = c,(0) 6 {0,1} 
the time-discrete dynamics is obtained by iterating the 
following update steps: 

(1) Choose a node i at random. (2) Calculate Si ac- 
cording to Eq. ©. (3) Update q according to Eq. JSJ). 

For t = and e = 1 random asynchronous update is 
recovered. For r > there is an explicit transmission 
delay from the output of node j to the input of node i. 
To be definite, at t — we assume that nodes have not 
flipped during the previous r time steps. 

Let us first explore the dynamics of a simple but non- 
trivial interaction network with N — 3 sites and non- 
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FIG. 2: Simulation runs of the three-node network in Fig.0 
(a) Random asynchronous update mode, r = 0, e = 1. (b) 
Filtering e = 0.01 but no delay r = 0. (c) Delay r = 100 
and no filtering, e = 1. A circle plotted at coordinates (t, i) 
indicates that state Xi of node i changes at time t. 
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FIG. 1: Left: Network with iV = 3 nodes and bias values 
ao = a,2 — and ai = 2. Right: Dynamics of the network. 
Transitions between configurations under asynchronous up- 
date are indicated by arrows. Under synchronous (parallel) 
update the system has one unique cyclic attractor only, con- 
sisting of the four configurations in the middle row. 

vanishing couplings woi — 1021 = — 1 and w\q — w%2 = 
+1, see Fig. ^ Note that under asynchronous update 
the sequence of states reached by the dynamics is not 
unique. The system may branch off to different config- 
urations depending on node update ordering. This is 
illustrated in Fig. |2{a) : Without delay (r = 0) and filter 
(e = 1) the dynamics is irregular, i.e. non-periodic. With 
filter only (r = 0, e = 0.01, Fig. EJb)), the dynamics is 
periodic at times, but also intervals of fast irregular flip- 
ping occur. Finally, in the presence of delay (r = 100, 
e = 1, Fig. GIc)) we obtain perfectly ordered dynam- 
ics with synchronization of flips. Nodes and 2 change 
states practically at the same (macro) time, followed by 
a longer pause until node 1 changes state, etc. With in- 
creasing delay time r the dynamics under asynchronous 
update approaches the dynamics under synchronous up- 
date (cf. Fig.^) when viewed on a coarse-grained (macro) 
time scale. 

Let us further quantify the difference between syn- 
chronous and asynchronous dynamics. First, a defini- 
tion of equivalence between the two dynamical modes 



has to be given. Let us start from the time series 
x(t) of configurations x = (xo, . . . , xn-i) produced by 
the asynchronous (random serial) update of the model 
and the respective time series y{u) produced by syn- 
chronous (parallel) update, using identical initial condi- 
tion y(0) = x(0). These time series live on different time 
scales, which we call the micro time scale of single site 
updates in the asynchronous case, and the macro time 
scale where each time step is an entire sweep of the sys- 
tem. Assume that at time t u the asynchronous system 
is in state x(t u ) = y(u). In order to follow the syn- 
chronous update it has to subsequently reach the state 
y(u +1) on a shortest path in phase space. Formally, 
let us require that there is a micro time t u+ i > t u such 
that x(t u +i) = y(u + 1) and each node flips at most 
once in the time interval [t u ,t u +i]. Once this is vio- 
lated we say that an error has occured at the partic- 
ular macro time step u. This error allows to define a 
numerical measure of discrepancy between asynchronous 
and synchronous dynamics. Starting from identical ini- 
tial conditions, the system is iterated in synchronous and 
asynchronous modes (here for utotai = 10 7 macro time 
steps). Whenever the resulting time series are no longer 
equivalent, an error counter is incremented and the sys- 
tem reset to initial condition. The total error E of the 
run is the number of errors divided by Mtotai- 

For the network in Fig. ^ and the initial condition 
Xi = Ci = for i = 1,2,3 the error E is exponentially 
suppressed with delay time r (Fig.|3|). The asynchronous 
dynamics with delay follows the attractor during a time 
span that increases exponentially with the given delay 
time. Note that there is only one possibility for the asyn- 
chronous dynamics to leave the attractor: When the sys- 
tem is in configuration (1,1,0) or (0,1,1), node 2 may 
change state such that the system goes to configuration 
(1,0,0) or (0,0, 1) respectively, whereas the correct next 
configuration on the attractor is (0,1,0). Consider the 
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FIG. 3: Discrepancy E between the asynchronous and syn- 
chronous update mode as a function of the delay time r and 
without filter (e = 1). The solid line is the theoretical predic- 
tion of the decay <x (2/3) r . 



case e = 1 where Cj = Xi for all i. Let us assume that the 
system is in configuration (1,1,1) and at time to node 
changes state, thereby generating configuration (0, 1, 1). 
This decreases the input sum si below zero such that for 
t = node would change state immediately in its next 
update. With explicit transmission delay r > 0, however, 
node 1 still "sees" the input sum Si = generated by the 
configuration (1, 1, 1) until time step to + r. If node 2 is 
chosen for update in this time window to + 1, ■ ■ ■ , to + r 
it changes state immediately and updates are performed 
in correct order. The opposite case, that node 2 does not 
receive an update in any of the r time steps, happens 
with probability (2/3) T , yielding the correct error decay 
of the simulation (Fig. [2J . 

Next we demonstrate that there are cases where also 
low-pass filtering, e <C 1, is needed for the asynchronous 
dynamics to follow the deterministic attractor. Consider 
a network of N = 5 nodes with bias values ao = a± = 
and ai = a-i = = 1. The only non-zero couplings are 
wio = Wn = ^31 = W42 = +1 and Woi = w 43 = — 1. 
Nodes and 1 form an oscillator, i.e. (xq, x\) iterate the 
sequence (0, 0), (1, 0), (1, 1), (0, 1). Nodes 2 and 3 simply 
"copy" the state of node 1 such that under synchronous 
update always x${t) — X2(t) — x\(t — 1). Consequently, 
under synchronous update the input sum of node 4 never 
changes because the positive contribution from node 2 
and the negative contribution from node 3 cancel out. 
Under asynchronous update, however, the input sum of 
node 4 may fluctuate because nodes 2 and 3 do not flip 
precisely at the same time. The effect of the low-pass 
filter e <C 1 is to suppress the spreading of such fluctua- 
tions on the micro time scale. The influence of the filter is 
seen in Fig.^J When r is kept constant, the error drops 
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FIG. 4: Discrepancy E between the asynchronous and syn- 
chronous update mode as a function of the filter parameter 
1/e for the network with N = 5 nodes described in the text. 
The delay parameter is chosen as r — (circles), t — 10 
(squares), and r = 1/e (diamonds). 



algebraically with decreasing e. An exponential decay 
E ~ exp(— a/e) is obtained when r oc 1/e (the filter can 
take full effect only in the presence of sufficient delay). 

Let us finally consider an example of a larger network 
with N = 16 nodes and L — 48 non-vanishing couplings 
(chosen randomly from the off-diagonal elements in the 
matrix (ife) and assigned values +1 or —1 with probabil- 
ity 1/2 each; biases are chosen as Oj = Y^j ?% /2). Sim- 
ulation runs under pure asynchronous update (r = 0, 
e = 1) typically yield dynamics as in Fig. Efa). The 
time series x(£) is non-periodic and non-reproducible, i.e. 
under different order of updates a different series is ob- 
tained. For the same initial condition, periodic dynam- 
ics is observed in the presence of sufficent transmission 
delay and filtering, Fig. 0Td). In this case, the system 
follows precisely the attractor of period 28 found under 
synchronous update. As seen in Fig. Etc), the error de- 
cays exponentially as a function of the delay time r. 

Let us now turn to the dangers of asynchronous up- 
date: There is a fraction of attractors observed under 
synchronous update that cannot be realized under asyn- 
chronous update. Synchronization cannot be sustained if 
the dynamics is separable. In the trivial case, separabil- 
ity means that the set of nodes can be divided into two 
subsets that do not interact with each other. Then there 
is no signal to synchronize one set of nodes with the other 
and they will go out of phase. In general, synchronization 
is impossible if the set of flips itself is separable. Con- 
sider, as the simplest example, a network of N = 2 nodes 
with the couplings wqi = wio = +1, biases ao = <X\ = 1 
and the initial condition (y (0), j/i(0)) = (0,1). Under 
synchronous update, the state alternates between vector 
(0,1) and (1,0). Under asynchronous update with de- 
lay time r, the transition of one node i from xi = to 
Xi = 1 causes the other node j to switch from Xj = 
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FIG. 5: Time series and errors obtained for the network 
of N = 16 nodes described in the text, (a) Asynchronous 
update mode with neither delay nor filtering, r = 0, e = 1. 
A vertical stroke at coordinates (t, i) indicates that node i 
flips at time t. At t « 1600 the system reaches a fixed point, 
(b) Same initial condition as in (a), but delay r = 50000 
and filtering e = 1/500. The system follows a limit cycle 
of 28 macro time steps, (c) Discrepancy (error E) between 
asynchronous and synchronous update mode as a function of 
the delay parameter for e = 50/r (circles), e = 40/t (squares), 
e = 25/r (diamonds). 



to Xj — 1 approximately r time steps later. The "on"- 
transitions only trigger subsequent "on" -transitions and, 
analogously, the "off " -transitions only trigger subsequent 
"off" -transitions. The dynamics can be divided into two 
distinct sets of events that do not influence each other. 
Consequently, synchronization between flips cannot be 
sustained, as illustrated in Fig. When the phase dif- 
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FIG. 6: Dynamics of the network with N = 2 mutually 
activating nodes. Rising edges and falling edges desynchronize 
over time. The two classes of edges form two uncoupled chains 
of events. 



ference reaches the value r, on- and off-transitions anni- 
hilate. Then the system leaves the attractor and reaches 
one of the fixed points with xq = x%. 



These observations have important implications for ro- 
bust topological motifs in asynchronous networks. First 
of all, the above example of a small excitatory loop can 
be quickly generalized to any larger loop with excitatory 
interactions, as well as to loops with an even number of 
inhibitory couplings, where in principle similar dynamics 
could occur. Higher order structures that fail to syn- 
chronize include competing modules, e.g. two oscillators 
(loops with odd number of inhibitory links) that interact 
with a common target. 



In conclusion we find that asynchronously updated net- 
works of autonomous dynamical nodes are able to ex- 
hibit a reproducible and quasi-deterministic dynamics 
under broad conditions if the nodes have transmission 
delay and low pass filtering as, e.g., observed in regu- 
latory genes. Timing requirements put constraints on 
the topology of the networks (e.g. suppression of certain 
loop motifs). With respect to biological gene regulation 
networks where indeed strong suppression of loop struc- 
tures is observed 0, 0] , one may thus speculate about 
a new constraint on topological motifs of gene regula- 
tion: The requirement for deterministic dynamics from 
asynchronous dynamical networks. 
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